############################################################
###########################Prediction code_ UNRB
############################################################
C:\Users\jwmille7\Dropbox\PC\Documents\R\win-library\4.1/00LOCK

####Set working directory (example:  setwd("K:/My Drive/EPA datasets"))
setwd("K:/")
setwd("Q:/")
setwd("C:/Users/jwmille7/Dropbox/EPA_WQ_Grant_Group_Files/Water Quality Data and Metadata/Final deliverables_ 07_2021/EPA datasets/")


library(lme4)

#Clear all variables
rm(list=ls())

set.seed=432
######Import models and data to recreate best models for wq parameters
######################################################
###Turbidity model
data.Tdu=readRDS("Prediction models/data.TDU_daily.rds")
wqdata.m=data.Tdu
fml.Tdu<-readRDS("Prediction models/fml.TDU.daily.rds")       # loads object "fml.As": list of best model formulas
Tdu.mod<-lmer(formula=fml.Tdu,data=data.Tdu)
summary(Tdu.mod)

# ###TN model
 data.TN=readRDS("Prediction models/data.TN_daily.rds")
 wqdata.m=data.TN
 fml.TN<-readRDS("Prediction models/fml.TN.daily.rds")       # loads object "fml.As": list of best model formulas
 TN.mod<-lmer(formula=fml.TN,data=data.TN)
 summary(TN.mod)
# 
# ###################################################################################
# ######Fix these three below#####
# 
# ########GEt these right in data files
# ###Fecal coliform model
 data.FC=readRDS("Prediction models/data.FC_daily.rds")
 wqdata.m=data.FC
 fml.FC<-readRDS("Prediction models/fml.FC.daily.rds")       # loads object "fml.As": list of best model formulas
 FC.mod<-lmer(formula=fml.FC,data=data.FC)
 summary(FC.mod)

# ###Macroinvertebrate model
 data.BI=readRDS("Prediction models/data.BI_yearly.rds")
 data.p=data.BI
 fml.BI<-readRDS("Prediction models/fml.BI.yearly.rds")       # loads object "fml.As": list of best model formulas
 BI.mod<-lmer(formula=fml.BI,data=data.BI)
 summary(BI.mod)


# ###TP model
data.TP=readRDS("Prediction models/data.TP_daily.rds")
 data.p=data.TP
 fml.TP<-readRDS("Prediction models/fml.TP.daily.rds")       # loads object "fml.As": list of best model formulas
 data.TP.tr=data.TP[complete.cases(data.TP),]
 names(data.TP)
 #TP.mod<-lmer(formula=fml.TP,data=data.TP)
 TP.mod<- lmer(data.TP[,1] ~ Lossofcc_30buffer
                       +Impervious_ws2
                       +log(TP_load_drainage_yr+1)
                       +log(X1WPrec+1) 
                       +X6MSPI
                       +Season #+Soil#+Year
                       + (1|SiteID)+(1|Basin2)+(1|Soil), data = data.TP,na.action=na.exclude,REML=F)
 summary(TP.mod)
# 
# ###Sp cond model
 data.SC=readRDS("Prediction models/data.SC_daily.rds")
# data.p=data.BI
  fml.SC<-readRDS("Prediction models/fml.SC.daily.rds")       # loads object "fml.As": list of best model formulas
  
  ###This is not working for some reason????
 # SC.mod<-lmer(formula=fml.SC,data=data.SC)
 
 SC.mod= lmer(data.SC[,1] ~ log(Lossofcc_30buffer+1)
               +Crop_ws
               #+log(Impervious_ws+1)
               +log(Impervious_ws5+1)
               +log(Impervious_wspre80.n+1)
               +log(Wetlands_30buffer+1)
               #+Wastewater_d
               #+log(TN_load_drainage_yr+1)
               +Majwwtp_num #+Lake5
               #+Minwwtp_num
               +log(X1WPrec+1)#+log(X1WPrec+1) 
               #+Impervious_ws:Lossofcc_30buffer
               #+Crop_ws:Lossofcc_30buffer
               #+Pasture_ws:Lossofcc_30buffer
               #+(log(X3DPrec+1)):Lossofcc_30buffer
               +X6MSPI
               #+Year
               +Season 
               + (1|SiteID)+(1|Basin2)+(1|Soil),data = data.SC,na.action=na.exclude,REML=F)
 summary(SC.mod)
 dat=model.frame(SC.mod)
 cor(dat[,c(1:9)])
 
 ###Extract all models and data
 t=Sys.Date()
write.csv(data.BI,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/BI_data",t,".csv",sep=""))
write.csv(data.Tdu,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/TDU_data",t,".csv",sep=""))
write.csv(data.FC,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/FC_data",t,".csv",sep=""))
write.csv(data.TN,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/TN_data",t,".csv",sep=""))
write.csv(data.TP,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/TP_data",t,".csv",sep=""))
write.csv(data.SC,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/SC_data",t,".csv",sep=""))
###Prediction datasets (must run model first in order to export these here)
write.csv(pred.sites,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/Predsites_normal_2",t,".csv",sep=""))
write.csv(pred.sites.m5,paste("EPA urban stream project-water quality/UNRB/EPA Datasets/Predsites_25pct red_cc_urban",t,".csv",sep=""))

 

#################################################
##############Making predictions for Upper Neuse sites using  models
##################################################
pred.sites= read.csv("Prediction models/NWALT_UNRB_newdata_2018-10-01.csv",header=TRUE,na.strings = "-9999")
colnames(pred.sites)[2]="SiteID"
#pred.sites=pred.sites[,c(2:28)] ##To get rid of first column...useless...
trial=pred.sites[pred.sites$SiteID=="A8778261",]

#######Adding reservoirs and WWTPs
pred.sites.wwtp= read.csv("Prediction models/UNRB_wwtp_lakes.csv",header=TRUE,na.strings = "-9999")
pred.sites$Lake=pred.sites.wwtp$Reservoir.5.km[match(pred.sites$SiteID,pred.sites.wwtp$SiteID)]
pred.sites$wwtpmaj.5=pred.sites.wwtp$within.1.km[match(pred.sites$SiteID,pred.sites.wwtp$SiteID)]

pred.sites$Lossofcc_30buffer_log= log(pred.sites$Lossofcc_30buffer+1)
pred.sites$Season="Summer"
pred.sites$Sample_Type="Full Scale"

hist(pred.sites$Impervious_ws5)
median(pred.sites$Impervious_ws5)
quantile(pred.sites$Impervious_ws5,.95)

#########GEt the correct TN load for each site for TN and TP models
library(reshape)
library(dplyr)

###Get TN, TP, and Discharge data sets
#################################################

wwtp.TN_yrmed = read.csv("Prediction models/WWTP_TN_by_year.csv",header=TRUE,na.strings = "-9999")
years=seq(1994,2016,by=1)
wwtp.TN_yrmed=wwtp.TN_yrmed[,c(2,4:26)];colnames(wwtp.TN_yrmed)[2:24]=years ###Get columns as numbers
TN.m=melt(wwtp.TN_yrmed);TN.m$variable=as.numeric(as.character(TN.m$variable)) ##Melt discharge
colnames(TN.m)=c("WWTPid","Year","TN_yrmed") 
TN.m=TN.m[TN.m$Year=="2012",]

wwtp.disch_yrmed = read.csv("Prediction models/WWTP_flow_by_year.csv",header=TRUE,na.strings = "-9999")
wwtp.disch_yrmed=wwtp.disch_yrmed[,c(2,4:26)];colnames(wwtp.disch_yrmed)[2:24]=years ###Get columns as numbers
flow.m=melt(wwtp.disch_yrmed);flow.m$variable=as.numeric(as.character(flow.m$variable)) ##Melt discharge
colnames(flow.m)=c("WWTPid","Year","Flow_yrmed") 

wwtp.TP_yrmed = read.csv("Prediction models/WWTP_TP_by_year.csv",header=TRUE,na.strings = "-9999")
wwtp.TP_yrmed=wwtp.TP_yrmed[,c(2,4:26)];colnames(wwtp.TP_yrmed)[2:24]=years ###Get columns as numbers
TP.m=melt(wwtp.TP_yrmed);TP.m$variable=as.numeric(as.character(TP.m$variable)) ##Melt discharge
colnames(TP.m)=c("WWTPid","Year","TP_yrmed") 
TP.m=TP.m[TP.m$Year=="2012",]


####Add yearly TN and TP loads for 2012 for all upstream Major and minor WWTPs
################################################################
pred.sites.wwtp$WWTP_Major1=TN.m$TN_yrmed[match(pred.sites.wwtp$WWTP_maj,TN.m$WWTPid)]
pred.sites.wwtp$WWTP_Major1[is.na(pred.sites.wwtp$WWTP_Major1)] <- 0 
pred.sites.wwtp$WWTP_Major1.p=TP.m$TP_yrmed[match(pred.sites.wwtp$WWTP_maj,TP.m$WWTPid)]
pred.sites.wwtp$WWTP_Major1.p[is.na(pred.sites.wwtp$WWTP_Major1.p)] <- 0 

pred.sites.wwtp$WWTP_Major2=TN.m$TN_yrmed[match(pred.sites.wwtp$Major,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Major2[is.na(pred.sites.wwtp$WWTP_Major2)] <- 0 
pred.sites.wwtp$WWTP_Major2.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Major,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Major2.p[is.na(pred.sites.wwtp$WWTP_Major2.p)] <- 0 

pred.sites.wwtp$WWTP_Minor1=TN.m$TN_yrmed[match(pred.sites.wwtp$Minor,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor1[is.na(pred.sites.wwtp$WWTP_Minor1)] <- 0 
pred.sites.wwtp$WWTP_Minor1.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Minor,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor1.p[is.na(pred.sites.wwtp$WWTP_Minor1.p)] <- 0 

pred.sites.wwtp$WWTP_Minor2=TN.m$TN_yrmed[match(pred.sites.wwtp$Minor.1,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor2[is.na(pred.sites.wwtp$WWTP_Minor2)] <- 0 
pred.sites.wwtp$WWTP_Minor2.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Minor.1,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor2.p[is.na(pred.sites.wwtp$WWTP_Minor2.p)] <- 0 

pred.sites.wwtp$WWTP_Minor3=TN.m$TN_yrmed[match(pred.sites.wwtp$Minor.2,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor3[is.na(pred.sites.wwtp$WWTP_Minor3)] <- 0 
pred.sites.wwtp$WWTP_Minor3.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Minor.2,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor3.p[is.na(pred.sites.wwtp$WWTP_Minor3.p)] <- 0 

pred.sites.wwtp$WWTP_Minor4=TN.m$TN_yrmed[match(pred.sites.wwtp$Minor.3,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor4[is.na(pred.sites.wwtp$WWTP_Minor4)] <- 0 
pred.sites.wwtp$WWTP_Minor4.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Minor.3,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor4.p[is.na(pred.sites.wwtp$WWTP_Minor4.p)] <- 0 

pred.sites.wwtp$WWTP_Minor5=TN.m$TN_yrmed[match(pred.sites.wwtp$Minor.4,TN.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor5[is.na(pred.sites.wwtp$WWTP_Minor5)] <- 0 
pred.sites.wwtp$WWTP_Minor5.p=TP.m$TP_yrmed[match(pred.sites.wwtp$Minor.4,TP.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor5.p[is.na(pred.sites.wwtp$WWTP_Minor5.p)] <- 0 


####Add yearly Discharge values for 2012 for all upstream Major and minor WWTPs
pred.sites.wwtp$WWTP_Major1.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$WWTP_maj,flow.m$WWTPid)]
pred.sites.wwtp$WWTP_Major1.f[is.na(pred.sites.wwtp$WWTP_Major1.f)] <- 0 

pred.sites.wwtp$WWTP_Major2.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Major,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Major2.f[is.na(pred.sites.wwtp$WWTP_Major2.f)] <- 0 

pred.sites.wwtp$WWTP_Minor1.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Minor,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor1.f[is.na(pred.sites.wwtp$WWTP_Minor1.f)] <- 0 

pred.sites.wwtp$WWTP_Minor2.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Minor.1,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor2.f[is.na(pred.sites.wwtp$WWTP_Minor2.f)] <- 0 

pred.sites.wwtp$WWTP_Minor3.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Minor.2,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor3.f[is.na(pred.sites.wwtp$WWTP_Minor3.f)] <- 0 

pred.sites.wwtp$WWTP_Minor4.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Minor.3,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor4.f[is.na(pred.sites.wwtp$WWTP_Minor4.f)] <- 0 

pred.sites.wwtp$WWTP_Minor5.f=flow.m$Flow_yrmed[match(pred.sites.wwtp$Minor.4,flow.m$WWTPid)] 
pred.sites.wwtp$WWTP_Minor5.f[is.na(pred.sites.wwtp$WWTP_Minor5.f)] <- 0 

#####################Join TN conc and flows to get TN loads/drainage
###Yearly medians TN
pred.sites.wwtp$TN_load_yr=(86.4*.0438)*(pred.sites.wwtp$WWTP_Major1*pred.sites.wwtp$WWTP_Major1.f 
                                    +pred.sites.wwtp$WWTP_Major2*pred.sites.wwtp$WWTP_Major2.f
                                    +pred.sites.wwtp$WWTP_Minor1*pred.sites.wwtp$WWTP_Minor1.f
                                    +pred.sites.wwtp$WWTP_Minor2*pred.sites.wwtp$WWTP_Minor2.f
                                    +pred.sites.wwtp$WWTP_Minor3*pred.sites.wwtp$WWTP_Minor3.f
                                    +pred.sites.wwtp$WWTP_Minor4*pred.sites.wwtp$WWTP_Minor4.f
                                    +pred.sites.wwtp$WWTP_Minor5*pred.sites.wwtp$WWTP_Minor5.f)

###Yearly medians TP
pred.sites.wwtp$TP_load_yr=(86.4*.0438)*(pred.sites.wwtp$WWTP_Major1.p*pred.sites.wwtp$WWTP_Major1.f 
                                         +pred.sites.wwtp$WWTP_Major2.p*pred.sites.wwtp$WWTP_Major2.f
                                         +pred.sites.wwtp$WWTP_Minor1.p*pred.sites.wwtp$WWTP_Minor1.f
                                         +pred.sites.wwtp$WWTP_Minor2.p*pred.sites.wwtp$WWTP_Minor2.f
                                         +pred.sites.wwtp$WWTP_Minor3.p*pred.sites.wwtp$WWTP_Minor3.f
                                         +pred.sites.wwtp$WWTP_Minor4.p*pred.sites.wwtp$WWTP_Minor4.f
                                         +pred.sites.wwtp$WWTP_Minor5.p*pred.sites.wwtp$WWTP_Minor5.f)

###Join to main dataset and divide by drainage area
pred.sites$TN_load_drainage_yr=pred.sites.wwtp$TN_load_yr/pred.sites$Drainage
pred.sites$TP_load_drainage_yr=pred.sites.wwtp$TP_load_yr/pred.sites$Drainage

####Add Majwwtp_num to predictors
pred.sites.wwtp$Maj_count1=ifelse(pred.sites.wwtp$WWTP_maj!=0,1,0)
pred.sites.wwtp$Maj_count2=ifelse(pred.sites.wwtp$Major!=0,1,0)
pred.sites$Majwwtp_num=pred.sites.wwtp$Maj_count1+pred.sites.wwtp$Maj_count2

####Change/add predictors of importance
pred.sites$Year=pred.sites$Year-2000
pred.sites$Minwwtp_num=0
pred.sites$Basin="No_basin"
pred.sites$Basin2="No_Basin"

#####To match predictor names and log values in newdata set
pred.sites$Impervious_ws80=pred.sites$Impervious_wspre80.n 
pred.sites$Impervious_ws80=log(pred.sites$Impervious_ws80+1)
pred.sites$Impervious_wspre80.n =log(pred.sites$Impervious_wspre80.n +1)
pred.sites$Impervious_ws=log(pred.sites$Impervious_ws+1)

##########################
####Getting 30 years of precipitation data for a site
##########################
###DomID for all prediction sites
domsitematch = read.csv("Prediction models/UNRB_all_domIDs.csv",header=TRUE,na.strings = "-9999")
pred.sites$domID=domsitematch$DomID[match(pred.sites$SiteID,domsitematch$SiteID)]

###Prec data 
domdata3 = read.csv("Prediction models/Jonathan_Piedmont_Precipitation_3day.csv",header=TRUE,na.strings = "-9999")
site.dom=as.character(pred.sites[1,"domID"]) 
prec=domdata3[,site.dom];prec.ts=cbind(domdata3[,1],prec) ##Gets only needed row
prec.ts=as.data.frame(prec.ts[c(28:length(prec.ts[,1])),]) ##To get 1980 to 2010 only
prec.ts$V1=prec.ts$V1-.0025
prec.ts$Year=floor(prec.ts[,1])
mat.year.3=split(prec.ts,prec.ts$Year)  ###To split up data by year for daily data
t=mat.year.3[[5]]  ###To see one year

###One week/ 7 day precipitation
domdata7 = read.csv("Prediction models/Jonathan_Piedmont_Precipitation_7day.csv",header=TRUE,na.strings = "-9999")
site.dom=as.character(pred.sites[1,"domID"]) 
prec=domdata7[,site.dom];prec.ts=cbind(domdata7[,1],prec) ##Gets only needed row
prec.ts=as.data.frame(prec.ts[c(28:length(prec.ts[,1])),]) ##To get 1980 to 2010 only
prec.ts$V1=prec.ts$V1-.0025
prec.ts$Year=floor(prec.ts[,1])
mat.year.7=split(prec.ts,prec.ts$Year)  ###To split up data by year for daily data
t=mat.year.7[[5]]  ###To see one year

###Two week/ 14 day precipitation
domdata14 = read.csv("Prediction models/Jonathan_Piedmont_Precipitation_14day.csv",header=TRUE,na.strings = "-9999")
site.dom=as.character(pred.sites[1,"domID"]) 
prec=domdata14[,site.dom];prec.ts=cbind(domdata14[,1],prec) ##Gets only needed row
prec.ts=as.data.frame(prec.ts[c(28:length(prec.ts[,1])),]) ##To get 1980 to 2010 only
prec.ts$V1=prec.ts$V1-.0025
prec.ts$Year=floor(prec.ts[,1])
mat.year.14=split(prec.ts,prec.ts$Year)  ###To split up data by year for daily data
t=mat.year.14[[5]]  ###To see one year

###############################################################
#####Function to get yearly sums of all 325 sites by year to see variation
prec.var=function(precdata){
        #precdata=domdata3
        ####Matrix to store data in
        cols=326-2+1; rows=31; 
        mat.3d.sums=as.data.frame(matrix(0,nrow=rows,ncol=cols))
        
        for (r in 2:326){
        #r=5
        prec=precdata[,r]; prec.ts=cbind(precdata[,1],prec) ##Gets only needed row
        prec.ts=as.data.frame(prec.ts[c(28:length(prec.ts[,1])),]) ##To get 1980 to 2010 only
        prec.ts$V1=prec.ts$V1-.0025
        prec.ts$Year=floor(prec.ts[,1])
        mat.year=split(prec.ts,prec.ts$Year)  ###To split up data by year for daily data
        
        ###################################
        #####Setting up matrix with 365 x 30 years of precipitation data
        
        lab=seq(1980,2010,1)
        cols=length(lab)+1; rows=365; 
        mat.3d=as.data.frame(matrix(0,nrow=rows,ncol=cols))
        mat.3d[,1]=seq(as.Date("1900/1/1"), as.Date("1900/12/31"), "day")
        colnames(mat.3d)=c("Day",lab)
        for (i in 1:31){
          dat=mat.year[[i]]
          mat.3d[,i+1]=dat[c(1:365),2]
             }
        
        ###To check summary statistics... TOO HIGH!!! 3 DAy sum.... 3x too high.
        means=as.data.frame(apply(mat.3d[,c(2:32)],2,mean,na.rm=T))  ###Change for wetlands
        sums=as.data.frame(apply(mat.3d[,c(2:32)],2,sum,na.rm=T))  ###Change for wetlands
        sums=sums/3000
        mat.3d.sums[,(r-1)]=sums
        }
  
  
  return(mat.3d.sums)
}

mat.3d.sums=prec.var(domdata3)
library(reshape)
t=melt(mat.3d.sums)
boxplot(t[,2],ylab="meters",main="Precipitation for all sites")

###################################
#####Setting up function to get matrix with 365 x 30 years of precipitation data
dailymat=function(mat.year){
lab=seq(1980,2010,1)
cols=length(lab)+1; rows=365; 
mat.3d=as.data.frame(matrix(0,nrow=rows,ncol=cols))
mat.3d[,1]=seq(as.Date("1900/1/1"), as.Date("1900/12/31"), "day")
colnames(mat.3d)=c("Day",lab)
for (i in 1:31){
  dat=mat.year[[i]]
  mat.3d[,i+1]=dat[c(1:365),2]
}
return(mat.3d)

}

mat.3d=dailymat(mat.year.3)  ##Get 3 day matrix for predictions
mat.7d=dailymat(mat.year.7)  ##Get 14 day matrix for predictions (IF needed)
mat.14d=dailymat(mat.year.14)  ##Get 14 day matrix for predictions (IF needed)


#############################################################
####################Drought data
#####Add something to account for North, south, central 
drought.n=  read.csv("Prediction models//Miller_Jonathan_ClimDiv3_Monthly.csv",header=TRUE,na.strings = "-9999")
drought.c=  read.csv("Prediction models//Miller_Jonathan_ClimDiv4_Monthly.csv",header=TRUE,na.strings = "-9999")

monthlymat=function(value,region){
  #######To get drought data
  z=melt(region,id = c("Year","Month"), measure.vars = value)
  z2=split(z,z$Year)
  

  ########Set up matrix 
  lab=seq(1980,2010,1)
  cols=length(lab)+1; rows=365; 
  mat.1spi=as.data.frame(matrix(0,nrow=rows,ncol=cols))
  mat.1spi[,1]=seq(as.Date("1900/1/1"), as.Date("1900/12/31"), "day")
  colnames(mat.1spi)=c("Day",lab) 
  mat.1spi$month=as.numeric(format(mat.1spi$Day, format = "%m"))
  for (l in 1:31){
    #l=3
    dat=z2[[l]]  ###Get monthly values for a given year
    months=mat.1spi[,c(1,33)]
    months$value=dat$value[match(months$month,dat$Month)]
    mat.1spi[,l+1]=months[,3]
    }
  return(mat.1spi[,c(1:32)])
}  ###Function to get drought data

mat.X1MSPI.n=monthlymat("X1MSPI",drought.n)
mat.X1MSPI.c=monthlymat("X1MSPI",drought.c)
mat.X6MSPI.n=monthlymat("X6MSPI",drought.n)
mat.X6MSPI.c=monthlymat("X6MSPI",drought.c)


#############################################
################Getting predictions
#############################################

###############################################
############Basic predictions for all models Dry/80th percentile weather and Management Scenario #1
###############################################
cols=25; rows=length(pred.sites[,2]);
pred.easy=as.data.frame(matrix(0,nrow=rows,ncol=cols))
pred.easy[,1]=pred.sites[,2]
lab.e1=c("BI","TDU","FC","TN","TP","SC");lab.e2=c("BI_wet","TDU_wet","FC_wet","TN_wet","TP_wet","SC_wet")
lab.m1=c("BI_m1","TDU_m1","FC_m1","TN_m1","TP_m1","SC_m1","BI_wet_m1","TDU_wet_m1","FC_wet_m1","TN_wet_m1","TP_wet_m1","SC_wet_m1")
colnames(pred.easy)=c("SITEID",lab.e1,lab.e2,lab.m1) 

######Set for average long-term rainfall and no short-term rainfall
per=.50
X3D=as.matrix(mat.3d[,c(2:32)]);X3Dper=quantile(X3D,per)/100  ###Percentile of 3day rain
X7D=as.matrix(mat.7d[,c(2:32)]);X7Dper=quantile(X7D,per)/100  ###Percentile of 3day rain
X14D=as.matrix(mat.14d[,c(2:32)]);X14Dper=quantile(X14D,per)/100  ###Percentile of 3day rain
X1MP=quantile(drought.n$X1MPrec,per,na.rm=T)*2.54/100
X3MP=quantile(drought.c$X3MPrec,per,na.rm=T)*2.54/100

pred.sites$X3MPrec=X3MP; pred.sites$X1MPrec=X1MP
pred.sites$X3MPrecsq=X3MP^2; pred.sites$X1MPrecsq=X1MP^2
pred.sites$X3DPrec=0;pred.sites$X1WPrec=0; pred.sites$X2WPrec=0
pred.sites$X1MSPI=0; pred.sites$X6MSPI=0; pred.sites$X12MSPI=0; pred.sites$X24MSPI=0
pred.sites$X3DPrec=0

#####Getting newdata for wet predictions 90th percentiles.... 
per=.90
pred.sites.w=pred.sites ###Start the newdata for wet conditions
X3D=as.matrix(mat.3d[,c(2:32)]);X3Dper=quantile(X3D,per)/100  ###Percentile of 3day rain
X7D=as.matrix(mat.7d[,c(2:32)]);X7Dper=quantile(X7D,per)/100  ###Percentile of 3day rain
X14D=as.matrix(mat.14d[,c(2:32)]);X14Dper=quantile(X14D,per)/100  ###Percentile of 3day rain
X1MP=quantile(drought.n$X1MPrec,per,na.rm=T)*2.54/100
X3MP=quantile(drought.n$X3MPrec,per,na.rm=T)*2.54/100
X1MSPI=quantile(drought.n$X1MSPI,per,na.rm=T)
X3MSPI=quantile(drought.n$X3MSPI,per,na.rm=T)
X6MSPI=quantile(drought.n$X6MSPI,per,na.rm=T)
X12MSPI=quantile(drought.n$X12MSPI,per,na.rm=T)
X24MSPI=quantile(drought.n$X24MSPI,per,na.rm=T)


####Wet current conditions
pred.sites.w$X3DPrec=X3Dper; pred.sites.w$X1WPrec=X7Dper; pred.sites.w$X2WPrec=X14Dper
pred.sites.w$X1MPrec=X1MP; pred.sites.w$X1MPrecsq=X1MP^2 
pred.sites.w$X3MPrec=X3MP; pred.sites.w$X3MPrecsq=X3MP^2 
pred.sites.w$X1MSPI= X1MSPI;pred.sites.w$X3MSPI= X3MSPI
pred.sites.w$X6MSPI= X6MSPI;pred.sites.w$X12MSPI= X12MSPI;pred.sites.w$X24MSPI= X24MSPI


#####Normalize drainage area in predictor datasets (mean conditions and wet conditions). 
hist(pred.sites$Drainage)
pred.sites$Drainage=pred.sites$Drainage-155.4   ##Center data (FC model)
pred.sites$Drainage=pred.sites$Drainage/134.3   ##Divide by standard deviation (FC model)
hist(pred.sites$Drainage)

hist(pred.sites.w$Drainage)
pred.sites.w$Drainage=pred.sites.w$Drainage-155.4   ##Center data (FC model)
pred.sites.w$Drainage=pred.sites.w$Drainage/134.3   ##Divide by standard deviation (FC model)
hist(pred.sites.w$Drainage)


######Management scenario #1 - Combination Canopy cover buffer +25%,impervious -25%
###Change prediction dataset based on management scenario for normal and wet conditions
summary(FC.mod)
pred.sites.m1=pred.sites

##Canopy loss 525
hist(pred.sites.m1$Lossofcc_30buffer)
hist(pred.sites.m1$Lossofcc_30buffer_log)
pred.sites.m1$Lossofcc_30buffer=0.75*pred.sites.m1$Lossofcc_30buffer
pred.sites.m1$Lossofcc_30buffer_log=log(pred.sites.m1$Lossofcc_30buffer+1)

#####Impervious impact 25%
hist(pred.sites.m1$Impervious_ws)
hist(pred.sites.m1$Impervious_wspre80.n)
pred.sites.m1$Impervious_ws=.75*pred.sites.m1$Impervious_ws
pred.sites.m1$Impervious_wspre80.n=.75*pred.sites.m1$Impervious_wspre80.n
pred.sites.m1$Impervious_wspre00.n=.75*pred.sites.m1$Impervious_wspre00.n
pred.sites.m1$Impervious_30buffer=.75*pred.sites.m1$Impervious_30buffer
pred.sites.m1$Impervious_ws5=.75*pred.sites.m1$Impervious_ws5


####Wet predictions Management scenario 1
pred.sites.w.m1=pred.sites.m1
pred.sites.w.m1$X3DPrec=X3Dper; pred.sites.w.m1$X1WPrec=X7Dper; pred.sites.w.m1$X2WPrec=X14Dper
pred.sites.w.m1$X1MPrec=X1MP; pred.sites.w.m1$X1MPrecsq=X1MP^2 
pred.sites.w.m1$X3MPrec=X3MP; pred.sites.w.m1$X3MPrecsq=X3MP^2 
pred.sites.w.m1$X1MSPI= X1MSPI;pred.sites.w.m1$X3MSPI= X3MSPI
pred.sites.w.m1$X6MSPI= X6MSPI;pred.sites.w.m1$X12MSPI= X12MSPI;pred.sites.w$X24MSPI= X24MSPI

##########################################
###########Loop for MEAN PREDICTIONS FOR ALL SCENARIOS
#############################################

ptm=proc.time()
for (s in 1:length(pred.sites[,1]) ){
  #s=787 ###length(pred.sites[,1])
  
  ##########Current conditions#########################
  summary(TN.mod)
  ###Dry predictions
  pred.bi=predict(BI.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  pred.tdu=predict(Tdu.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  pred.fc=predict(FC.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  pred.tn=predict(TN.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  pred.tp=predict(TP.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  pred.sc=predict(SC.mod,newdata=pred.sites[s,],allow.new.levels=TRUE)
  
  #####Wet predictions
  pred.bi.w=predict(BI.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  pred.tdu.w=predict(Tdu.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  pred.fc.w=predict(FC.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  pred.tn.w=predict(TN.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  pred.tp.w=predict(TP.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  pred.sc.w=predict(SC.mod,newdata=pred.sites.w[s,],allow.new.levels=TRUE)
  
  ###All current Predictions into matrix
  pred.easy[s,2]=pred.bi; pred.easy[s,8]=pred.bi.w
  pred.easy[s,3]=exp(pred.tdu)-1; pred.easy[s,9]=exp(pred.tdu.w)-1
  pred.easy[s,4]=exp(pred.fc)-1; pred.easy[s,10]=exp(pred.fc.w)-1
  pred.easy[s,5]=exp(pred.tn)-1; pred.easy[s,11]=exp(pred.tn.w)-1
  pred.easy[s,6]=exp(pred.tp)-1; pred.easy[s,12]=exp(pred.tp.w)-1
  pred.easy[s,7]=exp(pred.sc)-1; pred.easy[s,13]=exp(pred.sc.w)-1
  
  ###########MAnagement Scenario 1 Combination all IC and Canopy Cover improvements#########################
  #summary(SC.mod)
  library(merTools)
  predictInterval(BI.mod,newdata=pred.sites[s,])
  predictInterval(BI.mod,newdata=pred.sites[s,],which="fixed")
  predictInterval(BI.mod,newdata=pred.sites[s,],which="all")
  predictInterval(FC.mod,newdata=pred.sites[s,])
  predictInterval(FC.mod,newdata=pred.sites[s,],which="fixed")
  predictInterval(FC.mod,newdata=pred.sites[s,],which="all")
  ###Mean predictions
  pred.bi=predict(BI.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  pred.tdu=predict(Tdu.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  pred.fc=predict(FC.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  pred.tn=predict(TN.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  pred.tp=predict(TP.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  pred.sc=predict(SC.mod,newdata=pred.sites.m1[s,],allow.new.levels=TRUE)
  
  #####Wet predictions
  pred.bi.w=predict(BI.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  pred.tdu.w=predict(Tdu.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  pred.fc.w=predict(FC.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  pred.tn.w=predict(TN.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  pred.tp.w=predict(TP.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  pred.sc.w=predict(SC.mod,newdata=pred.sites.w.m1[s,],allow.new.levels=TRUE)
  
  ###All current Predictions into matrix
  pred.easy[s,14]=pred.bi; pred.easy[s,20]=pred.bi.w
  pred.easy[s,15]=exp(pred.tdu)-1; pred.easy[s,21]=exp(pred.tdu.w)-1
  pred.easy[s,16]=exp(pred.fc)-1; pred.easy[s,22]=exp(pred.fc.w)-1
  pred.easy[s,17]=exp(pred.tn)-1; pred.easy[s,23]=exp(pred.tn.w)-1
  pred.easy[s,18]=exp(pred.tp)-1; pred.easy[s,24]=exp(pred.tp.w)-1
  pred.easy[s,19]=exp(pred.sc)-1; pred.easy[s,25]=exp(pred.sc.w)-1
  
  }
(proc.time()[3]-ptm[3])/60 #finish timing run

##Export file where you want it - Pred easy for Management Scenario #1 (25% cc and IC, 0% wwtp)
t=Sys.Date()
write.csv(pred.easy,paste("Prediction models/ManagementScenario1 Predictions_",t,".csv",sep=""))


######Management scenario #4 - decrease Random effects in sites and basins by 25%

####Setting up random initialization of random effects for each model
set.seed=444
##Get BI random effects
BI.re=as.data.frame(VarCorr(BI.mod)) 
BI.s.sd=BI.re[1,5]
BI.b.sd=BI.re[2,5]  ###Basin st. dev

summary(BI.mod)

####Length of prediction sites and percent reduction of random effects
s=length(pred.sites.m4[,1]);

#####Set the percent for the random effect reduction 1 = no reduction. 0 = 100% reduction
per=.75; ####25% reduction in positive random effects (those that degrade WQ)
#per=1

#BI model
m=0;sd=BI.s.sd;sd.b=BI.b.sd 

#red.BI=vector(,s); ####To get only one vector
cols=100; rows=s;
mat.BI=as.data.frame(matrix(0,nrow=rows,ncol=cols))

for (c in 1:cols){
  red.BI=vector(,s); red.BI.s=vector(,s);red.BI.b=vector(,s) ####To get vectors
  for (i in 1:s){
    re.s=rnorm(1,mean=m,sd=sd)
    re.b=rnorm(1,mean=m,sd=sd.b)
    
    if (re.s >0){
      red.BI.s[i]=re.s*per
    } else {
      red.BI.s[i]=re.s
    }
    if (re.b >0){
      red.BI.b[i]=re.b*per
    } else {
      red.BI.b[i]=re.b
    }
    red.BI=red.BI.s+red.BI.b
    mat.BI[,c]=red.BI
  }
}
mat.BI.sum=apply(mat.BI,2,mean)
mat.BI.mean=sum(mat.BI.sum)/100

##Get FC random effects
FC.re=as.data.frame(VarCorr(FC.mod)) 
FC.s.sd=FC.re[1,5]  ###Site st. dev
FC.b.sd=FC.re[2,5]  ###Basin st. dev

#FC model
m=0;sd.s=FC.s.sd;sd.b=FC.b.sd

cols=100; rows=s;
mat.FC=as.data.frame(matrix(0,nrow=rows,ncol=cols))

for (c in 1:cols){
  red.FC=vector(,s);red.FC.s=vector(,s);red.FC.b=vector(,s)
for (i in 1:s){
  #i=1
  re.s=rnorm(1,mean=m,sd=sd.s)
  re.b=rnorm(1,mean=m,sd=sd.b)
  if (re.s >0){
    red.FC.s[i]=re.s*per
  } else {
    red.FC.s[i]=re.s
  }
  if (re.b >0){
    red.FC.b[i]=re.b*per
  } else {
    red.FC.b[i]=re.b
  }
red.FC=red.FC.s+red.FC.b
}
  mat.FC[,c]=red.FC
  
  }
mat.FC.sum=apply(mat.FC,2,mean)
mat.FC.mean=sum(mat.FC.sum)/100

##Get TDU random effects
##Get TDU random effects
TDU.re=as.data.frame(VarCorr(Tdu.mod)) 
TDU.s.sd=TDU.re[1,5]  ###Site st. dev
TDU.b.sd=TDU.re[2,5]  ###Basin st. dev

cols=100; rows=s;
mat.TDU=as.data.frame(matrix(0,nrow=rows,ncol=cols))

#TDU model
m=0;sd.s=TDU.s.sd;sd.b=TDU.b.sd 

for (c in 1:cols){
  red.TDU=vector(,s);red.TDU.s=vector(,s);red.TDU.b=vector(,s)
  for (i in 1:s){
    #i=1
    re.s=rnorm(1,mean=m,sd=sd.s)
    re.b=rnorm(1,mean=m,sd=sd.b)
    if (re.s >0){
      red.TDU.s[i]=re.s*per
    } else {
      red.TDU.s[i]=re.s
    }
    if (re.b >0){
      red.TDU.b[i]=re.b*per
    } else {
      red.TDU.b[i]=re.b
    }
    red.TDU=red.TDU.s+red.TDU.b
  }
  mat.TDU[,c]=red.TDU
  
}

mat.TDU.sum=apply(mat.TDU,2,mean)
mat.TDU.mean=sum(mat.TDU.sum)/100
sum=sum(red.TDU);effect.TDU=sum/s

##Get TN random effects
TN.re=as.data.frame(VarCorr(TN.mod)) 
TN.s.sd=TN.re[1,5]  ###Site st. dev
TN.b.sd=TN.re[2,5]  ###Basin st. dev

#TN model
cols=100; rows=s;
mat.TN=as.data.frame(matrix(0,nrow=rows,ncol=cols))

#TN model
m=0;sd.s=TN.s.sd;sd.b=TN.b.sd 

for (c in 1:cols){
  red.TN=vector(,s);red.TN.s=vector(,s);red.TN.b=vector(,s)
  for (i in 1:s){
    #i=1
    re.s=rnorm(1,mean=m,sd=sd.s)
    re.b=rnorm(1,mean=m,sd=sd.b)
    if (re.s >0){
      red.TN.s[i]=re.s*per
    } else {
      red.TN.s[i]=re.s
    }
    if (re.b >0){
      red.TN.b[i]=re.b*per
    } else {
      red.TN.b[i]=re.b
    }
    red.TN=red.TN.s+red.TN.b
  }
  mat.TN[,c]=red.TN
  
}

mat.TN.sum=apply(mat.TN,2,mean)
mat.TN.mean=sum(mat.TN.sum)/100

sum=sum(red.TN);effect.TN=sum/s

##Get TP random effects
TP.re=as.data.frame(VarCorr(TP.mod)) 
TP.s.sd=TP.re[1,5]  ###Site st. dev
TP.b.sd=TP.re[2,5]  ###Basin st. dev

cols=100; rows=s;
mat.TP=as.data.frame(matrix(0,nrow=rows,ncol=cols))

#TP model
m=0;sd.s=TP.s.sd;sd.b=TP.b.sd 

for (c in 1:cols){
  red.TP=vector(,s);red.TP.s=vector(,s);red.TP.b=vector(,s)
  for (i in 1:s){
    #i=1
    re.s=rnorm(1,mean=m,sd=sd.s)
    re.b=rnorm(1,mean=m,sd=sd.b)
    if (re.s >0){
      red.TP.s[i]=re.s*per
    } else {
      red.TP.s[i]=re.s
    }
    if (re.b >0){
      red.TP.b[i]=re.b*per
    } else {
      red.TP.b[i]=re.b
    }
    red.TP=red.TP.s+red.TP.b
  }
  mat.TP[,c]=red.TP
  
}

mat.TP.sum=apply(mat.TP,2,mean)
mat.TP.mean=sum(mat.TP.sum)/100

sum=sum(red.TP);effect.TP=sum/s

##Get SC random effects
SC.re=as.data.frame(VarCorr(SC.mod)) 
SC.s.sd=SC.re[1,5]  ###Site st. dev
SC.b.sd=SC.re[2,5]  ###Basin st. dev

cols=100; rows=s;
mat.SC=as.data.frame(matrix(0,nrow=rows,ncol=cols))


#SC model
m=0;sd.s=SC.s.sd;sd.b=SC.b.sd 

for (c in 1:cols){
  red.SC=vector(,s);red.SC.s=vector(,s);red.SC.b=vector(,s)
  for (i in 1:s){
    #i=1
    re.s=rnorm(1,mean=m,sd=sd.s)
    re.b=rnorm(1,mean=m,sd=sd.b)
    if (re.s >0){
      red.SC.s[i]=re.s*per
    } else {
      red.SC.s[i]=re.s
    }
    if (re.b >0){
      red.SC.b[i]=re.b*per
    } else {
      red.SC.b[i]=re.b
    }
    red.SC=red.SC.s+red.SC.b
  }
  mat.SC[,c]=red.SC
  
}

mat.SC.sum=apply(mat.SC,2,mean)
mat.SC.mean=sum(mat.SC.sum)/100

 
#####mat.WQ.mean is the random effect adjustment needed to calculate MS2 from normal conditions 
####(taking into account log transformation except for BI)
####Based on 142,300 random draws

mat.BI.sum=apply(mat.BI,2,mean)
mat.BI.mean=sum(mat.BI.sum)/100

mat.FC.sum=apply(mat.FC,2,mean)
mat.FC.mean=sum(mat.FC.sum)/100

mat.TDU.sum=apply(mat.TDU,2,mean)
mat.TDU.mean=sum(mat.TDU.sum)/100

mat.TN.sum=apply(mat.TN,2,mean)
mat.TN.mean=sum(mat.TN.sum)/100

mat.TP.sum=apply(mat.TP,2,mean)
mat.TP.mean=sum(mat.TP.sum)/100

mat.SC.sum=apply(mat.SC,2,mean)
mat.SC.mean=sum(mat.SC.sum)/100
